In order of occurence:

• rgdal (reading in shapefiles)
• fields (color scheme)
• geosphere (calculating distance between points)
• tidyverse (mostly for ggplot)
• maps (for clean basemap)

• Any data with that contains geographic coordinates (GPS)
• Can be tables, shapefiles, and even photographs
• Example: looking at maps to see the most recent earthquake

• Locations are usually given by two distinct numbers

• To designate direction: sometimes positive or negative numbers, sometimes positive numbers with N/S or E/W next to them!

- To designate location: decimal degrees (DD), degrees minutes seconds (DMS), degrees minutes.minutes (DM.M), and sometimes even meters (what?!)

So… geospatial data can be tricky, and a bit sticky. But once you're used to it, it can be oh so satisfying!

• Tracking animals: Where do they go? How do they interact with their environment? (today)

• Tracking flight patterns: Can an airline add more flights to their schedule?
• Predicting the next earthquake from previous earthquake locations
• and so much more!

You're a marine biologist, studying fish movement patterns.

Research questions

• Are bottom-feeding fish spending time near a wastewater outfall?
• How close are these fish to the outfall pipe?

Implications:

• Fish could be exposed to harmful contaminants
• Humans that eat these fish could also be exposed

In your toolbox:

• shapefile with wastewater outfall location
• shapefile with water depths (bathymetry)
• shapefile with US states (land masses)
• tagged fish locations

Wait, what's a shapefile?

• A set of geospatial files (.shp, .shx, .dbf, .prj) that can map themselves
• The set must stay together, or the file breaks, but you only open the .shp

Step 1. Load required packages.

library(rgdal)

Step 2. Load shapefiles.

outfall <- readOGR("sample-data/shapefiles/outfall-reproject.shp",
                   verbose=FALSE)
depths <- readOGR("sample-data/shapefiles/all-depths.shp", verbose=FALSE)
states <- readOGR("sample-data/shapefiles/us-states.shp", verbose=FALSE)

For all the layers to plot nicely on top of each other, they need to all be in the same coordinate system.

Step 4. Convert coordinate systems.

states <- spTransform(states, proj4string(outfall))

Check that the projection has properly converted

proj4string(states)

## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Now, we can quickly plot the data and see how it looks so far.

plot(depths, col="blue", lty=2)
plot(outfall, col="black", lwd=2, add=TRUE)
plot(states, col="gray", add=TRUE)

Step 5: Load in fish and receiver GPS points.

receiver_locations <- read.csv("sample-data/csv/receiverlocs.csv")
fish_locations <- read.csv("sample-data/csv/location_data.csv")

head(receiver_locations, n=2)

##         Y         X Station   VR2W
## 1 33.5927 -118.0067  OCSD_1 129778
## 2 33.5909 -118.0018  OCSD_2 129789

head(fish_locations, n=2)

##   X               times       LON      LAT
## 1 1 2017-06-16 16:16:00 -117.9986 33.58484
## 2 2 2017-06-16 16:17:00 -117.9982 33.58450

For this plot, I want to be kind of fancy.
I want the fish locations to be color-coded by time.

# First, make sure date/time are in POSIX format
fish_locations$times <- as.POSIXct(fish_locations$times, tz="UTC",
                                   format="%Y-%m-%d %H:%M:%S")

# Then, convert to numeric so we can add a color scheme
fish_locations$numeric_time <- as.numeric(fish_locations$times)

Load in color scheme package.

library(fields)

Set plot boundaries (trial and error)

xlimits<-c(-118.076750, -117.953354) # only interested in this range of longitudes
ylimits<-c(33.558575, 33.637742) # and this range of latitudes

In depths shapefile, the "field_3" gives the depth in meters
Where -30 == 30 meters depth

library(geosphere)

# Grab random sample (for processing speed)
random_sample <- fish_locations[sample(nrow(fish_locations), 1000), ]

# Figure out the shapefile attribute that is closest to the fish gps point
distance_depths <- dist2Line(random_sample[,c("LON", "LAT")], depths)
distance_depths <- as.data.frame(distance_depths)

# Add depth column to random_sample that gives us the depth value
random_sample$depth <- depths$field_3[distance_depths$ID]

Plot it!

distance_table <- as.data.frame(table(random_sample$depth))
barplot(height=distance_table$Freq, names.arg=rev(distance_table$Var1),
        col="darkorange4", main="", xlab="depth (m)", ylab="frequency")

Calculate the distance between fish locations and the nearest part of the outfall pipe.

# Dist2Line to get distances
outfall_distance <- dist2Line(random_sample[,c("LON", "LAT")], outfall)

# Convert to data frame
outfall_distance <- as.data.frame(outfall_distance)

# Add distance (m) to the random sample data frame
random_sample$outfall_distance <- outfall_distance$distance

How do our data compare?

# Round fish locations into 250 m groups
random_sample$outfall_distance <- floor(random_sample$outfall_distance/250)*250

# Generate expected and observed plots, but don't plot them yet
expected <- hist(expected_hist, breaks=12, plot=FALSE)
obs <- hist(random_sample$outfall_distance, breaks=12, plot=FALSE)

So, from this, it looks like this individual isn't regularly using the outfall pipe.

Let's add some fish data and re-plot it.

# Add our fish data
layers <- worldplot + 
  geom_point(fish_locations, mapping = aes(x=LON, y=LAT), 
             color="darkorange", alpha=0.5)

# Create theme so background is blue like the ocean
new_theme <- theme(panel.background = 
                     element_rect(fill = "slategray1", colour = NA), 
        panel.grid.major = element_line(colour = NA),
        panel.grid.minor = element_line(colour = NA))

# Plot it
# layers + new_theme

• ggplot is a lot like matplotlib in python
• Lots of little steps to create a nice plot
• But it gets pretty tricky when you're trying to add your own shapefiles or trying to have a very fine-scale spatial resolution

Mapping in R can be difficult, but geospatial data can answer some very interesting questions.

Once you get the hang of it, mapping in R can be very rewarding!